close all;
load underneath gr well
load setup setup
gr.mesh;
axis vis3d equal tight;
% camlight;
% camlight(90,0);
lighting phong


%%  plot aquarium
glass.thickness     = 0.5;  % [cm]
glass.alpha         = 0.05;
glass.color         = [0.1 0.1 0.8];

% Bottom
plotcube([-max(gr.xGr) -max(gr.yGr) -glass.thickness], [max(gr.xGr) max(gr.yGr) min(gr.zGr)], glass.alpha, glass.color);    % (startcoordinates,size of cube,glassprop)
% Back
plotcube([-max(gr.xGr) glass.thickness abs(min(gr.zGr))], [max(gr.xGr) max(gr.yGr) min(gr.zGr)], glass.alpha, glass.color);
% RHS
plotcube([glass.thickness -max(gr.yGr) abs(min(gr.zGr))], [max(gr.xGr) max(gr.yGr) min(gr.zGr)], glass.alpha, glass.color);

% injectionbarrier
% horizontal
plotcube([-diff(setup(1).xhInjB) -diff(setup(1).yhInjB) abs(diff(setup(1).zhInjB))], [max(setup(1).xhInjB) max(setup(1).yhInjB) min(setup(1).zhInjB)], glass.alpha, glass.color)
% vertical
plotcube([-diff(setup(1).xvInjB) -diff(setup(1).yvInjB) abs(diff(setup(1).zvInjB))], [max(setup(1).xvInjB) max(setup(1).yvInjB) min(setup(1).zvInjB)], glass.alpha, glass.color);


%%  plot WEL
% [iyWEL  ixWEL   izWEL]  = ind2sub(size(IBOUND),find(IBOUND == 4));
well.plot3D('LineWidth',5);

%%  plot RECH
try
RECHactive  = sum(abs(RECH),3);
[iyRECH ixRECH izRECH]  = ind2sub(size(IBOUND),find(0 < RECHactive | RECHactive < 0));

RECHplot    = scatter3(gr.xm(ixRECH),gr.ym(iyRECH),gr.zm(izRECH));

catch
    if ~exist('RECH')
        mf_build;
        plottop;
    else
        error('')
    end
end

%%  plot HFB
ixHFB   = HFB(1,5);
iyHFB   = gr.yGr;
izHFB   = 1:(max(HFB(:,4))+1);

p1  = [gr.xGr(ixHFB) min(gr.yGr) max(gr.zGr(izHFB))];
p2  = [gr.xGr(ixHFB) max(gr.yGr) max(gr.zGr(izHFB))];
p3  = [gr.xGr(ixHFB) max(gr.yGr) min(gr.zGr(izHFB))];
p4  = [gr.xGr(ixHFB) min(gr.yGr) min(gr.zGr(izHFB))];

x   = [p1(1) p2(1) p3(1) p4(1)];
y   = [p1(2) p2(2) p3(2) p4(2)];
z   = [p1(3) p2(3) p3(3) p4(3)];

HFBplot     = fill3(x, y, z, 1);
alpha(HFBplot,0.5);


%% stub
% surf(peaks)
% axis vis3d
% for i=1:720
% 	camroll(.2)
% 	drawnow
% end

